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We consider the dynamics of Bianchi A scalar field models which undergo infla- 
tion. The main question is under which conditions does inflation come to an end and 
\C is succeeded by a decelerated epoch. This so-called "graceful exit" from inflation is 

'"^ an important ingredient in the standard model of cosmology, but is, at this stage, 

^_^ only understood for restricted classes of solutions. We present new results obtained 

O by a combination of analytical and numerical techniques. 

bJO 1 Introduction 

T— I Studies of cosmological solutions of Einstein's field equations have led to astonishing 

1^ insights about the history and fate of our own universe. In particular, the standard 

model of cosmology [24] has been remarkably successful, for example in describing the 
C<^ distribution of stars and galaxies over time, the primordial nucleosynthesis and the 

'^ cosmic microwave background (CMB). The standard model is based on the cosmological 

"^ principle. The idea is that there should be no preferred place nor direction in the universe 

Z^ and hence it should be spatially homogeneous and isotropic ~ at least on average over 

T-H large scales. In Einstein's theory of gravity, this idea leads to the FriedTnann- Roberts on- 

jj^ Walker (FRW) models. However, when FRW models are coupled to "normal" matter 

fields only and the results are compared with observations, several problems arise, for 
example the flatness problem and the horizon problem. The idea of inflation, which was 
Cd introduced for the first time in [14], addresses these problems. Inflation is a short phase 

of rapid expansion just after the big bang which is driven by a hypothetic matter field 
called inflaton. Indeed, inflation has become an important cornerstone of the standard 
model. An introduction to this idea from a rather mathematical point of view can be 
found in [28]. However, we also point the reader to the rather critical discussion in [8]. It 
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is generally believed that the inflationary epoch must have stopped after approximately 
75 — 100 e-folds and was then succeeded by an era of decelerated expansion where 
"normal" baryonic matter, radiation, neutrinos and dark matter became dominant. The 
transition from inflation to this decelerated stage is often called graceful exit. If a graceful 
exit from inflation exists then one also speaks of transient inflation as opposed to eternal 
inflation which lasts forever. In this paper, we focus on such transitions from inflation 
to a decelerated regime and hence the graceful exit problem. 

The standard model is a drastic simplification since our universe is clearly not pre- 
cisely spatially homogeneous and isotropic. Normally the idea is that within the light 
cone of an arbitrary observer, "initial" inhomogeneities and anisotropics, which are 
present just after the big bang, die out rapidly during inflation. Then, when infla- 
tion stops and the universe becomes decelerated, these should stay small in comparison 
to other measurable quantities. This kind of argument is often used as a justification 
why the standard model is expected to be a good description of the universe. How- 
ever, this is far from being clear. On the one hand, it is not clear how the averaging 
process which is necessary for describing the (presumably) approximately homogeneous 
and isotropic universe by exactly homogeneous and isotropic models is defined. In fact, 
it is possible that nonlinear effects give rise to back-reaction phenomena which are not 
taken into account by the standard model, see e.g. [10, 9, 18, 13]. On the other hand, 
the stability of Friedmann-Robertson- Walker solutions within the space of general so- 
lutions of Einstein's equation is not understood. Only in the case of certain eternally 
inflationary models, the problem of stability has been solved [31, 32]. It may therefore 
be the case that the standard model does not describe the whole history of our universe 
reliably. While there is relatively firm support for the idea that inhomogeneities and 
anisotropics indeed decay during inflation, it is particularly unclear what happens when 
inflation comes to an end. In order to address such fundamental questions, we need 
in principle to give up homogeneity and isotropy completely and study the dynamics 
of general classes of models. In practice it is, however, impossible to investigate com- 
pletely general solutions. In this paper here, we compromise by keeping the assumption 
of spatial homogeneity but give up isotropy. In this setting, the field equations imply 
"only" ordinary differential equations (plus algebraic constraints), which are relatively 
tractable. 

The simplest models for accelerated expansion are obtained by introducing a positive 
cosmological constant to Einstein's field equations^ and otherwise restrict to "normal" 
matter fields, i.e., matter fields which satisfy the dominant and strong energy conditions. 
A theorem by Wald [35] (which applies to the spatially homogeneous case) suggests that 
inflation is eternal for such models and there is therefore no graceful exit. Hence models 
with a cosmological constant cannot describe the transition from inflation to "ordinary" 
matter-dominated dynamics. In this paper, we study minimally coupled scalar field 
models. For such models, as we will describe in more detail later, we are free to choose the 
scalar field potential y as a smooth function of the scalar field (j). There is lots of freedom 

^We assume the signature ( — h++) for the spacetime metric. Our sign convention is that the de-Sitter 
spacetime is a solution of Einstein's vacuum equations with a positive cosmological constant. 



for this choice, but physical considerations give rise to some restrictions. The general 
idea is that, during the evolution of the universe, the scalar field rolls down the potential 
while experiencing a friction force whose magnitude is determined by the geometry of 
the spacetime. This suggests that the choice of the shape of the potential has important 
consequences for the evolution. Later on, we shall describe the current knowledge about 
the dynamics of Bianchi scalar field models. In turns out that most of the works on 
acceleratedly expanding scalar field models in the literature are targeted to situations 
where inflation does not stop and hence there is no graceful exit. This is so because 
the main interest of many authors lies in another fundamental outstanding problem of 
mathematical cosmology: the cosmic no-hair conjecture, which was formulated first in 
[12, 16], see also [6]. This conjecture states that generic expanding solutions of Einstein's 
field equations with reasonable matter fields and with a positive cosmological constant 
or a suitable scalar field approach the de-Sitter solution asymptotically. Nevertheless, 
even though these results exclude the actual case of our interest for this paper, namely 
inflation with a graceful exit, they give us valuable insights about the graceful exit 
problem. While we are aware of only a few works where the emphasis lies on the graceful 
exit (among those are [25, 4, 7], see also a simple example in [24]), most of which are 
restricted to the homogeneous and isotropic case, we consider this problem now in the 
anisotropic case. We choose a specific scalar field potential here (as discussed later) which 
to our knowledge has not been considered in this context before. We point out that in 
this paper, we are mainly interested in general qualitative properties of cosmological 
models as opposed to quantitative astrophysically relevant results. In order to make our 
discussion here as simple as possible and henceforth to be able to study our questions in 
a "clean" setting, we restrict to cosmological models where only one minimally coupled 
scalar field and no other matter fields are present. 

For our studies, the following strategy is adopted. We consider solutions of the 
Einstein-scalar field system from initial data in the infiationary regime. We then study 
the evolution of such models by a mixture of analytical and numerical techniques in order 
to determine under which conditions inflation only lasts a flnite time and hence there 
is a graceful exit. In turns out that for our choice of the potential and under certain 
further conditions, graceful exits from inflation occur generically due to the existence of 
certain future attractor solutions in the decelerated regime. 

The paper is organized as follows. In Section 2 we discuss the necessary background 
for Bianchi A scalar fleld models and derive the full set of equations implied by Ein- 
stein's equations and the scalar fleld equations. In Section 3 we describe our numerical 
scheme of choosing initial data and how to deal with numerical constraint violations. 
Section 4 is the main part of the paper. In Section 4.1, we summarize the main strategy. 
Basic consequences of our assumptions for the dynamics of the models are discussed in 
Section 4.2. Section 4.3 is devoted to a summary of the known results for various classes 
of potentials. Since all of these results exclude the case of inflation with a graceful exit, 
we find a potential which is not covered by these results. We then discuss the basic 
properties of this potential. It gives rise to three different cases. In one case, where the 
potential is strictly monotonically decreasing, we expect that the scalar field approaches 



oo asymptotically and a graceful exit is possible. Hence, we analyze this case in detail in 
Section 4.4 first. If the potential is not monotonic, but has rather a local minimum at a 
finite value of (p, it depends on the initial conditions whether the solution has a graceful 
exit. We study this case in Section 4.5. Finally, we close the paper with conclusions in 
Section 5. 

2 Background 

2.1 Self-gravitating minimally coupled scalar field models 

We consider globally hyperbolic, time-oriented oriented 4-dimensional smooth Lorentzian 
manifolds {M,g) where^ g = g^i, is a smooth Lorentzian metric. We are concerned 
with solutions of Einstein's field equation with vanishing cosmological constant (in ge- 
ometrized units c = SttG = 1 for the speed of light c and the gravitational constant 
G), 

where G^y is the Einstein tensor of g^y and T^i, is the energy momentum tensor of the 
matter fields. As the matter field, we choose a minimally coupled scalar field (p whose 
energy momentum tensor is 

V = V^'/'V^t/. - g^y (^Va^Vf^cj) + V{cj))\ . (2) 

The remaining freedom is the choice of the scalar field potential V{(j)) which is a smooth 
function of the scalar field (j). If this potential is strictly non-negative, as we assume 
throughout this paper, it follows that the dominant and hence weak energy conditions 
are satisfied. The strong energy condition, however, may be violated. Indeed, the 
violation of the strong energy condition is usually considered as the essential driving 
force of accelerated expansion and hence infiation. The equations of motion for the 
scalar field, 

V^V.,-^.0, (3) 

are derived from V^T^i^ = 0. Our particular choice of the function V{(j)) is discussed 
later. Notice that since no further matter fields are present, our models could also be 
described as self- gravitating scalar field models. 

2.2 Bianchi A spacetimes 

In this paper, we assume a particular class of spatially homogeneous models: Bianchi 
models. The basic assumption for Bianchi models is that the isometry group has a 
3-dimensional Lie subgroup which acts simply transitively by isometries on spacelike 

■^Our conventions for writing tensor fields are as follows. We either write the symbol of the tensor 
field without indices, e.g., g, or we use the abstract index notation, e.g., g^,,, with indices ^, v,. . . . Greek 
indices are therefore abstract indices and hence do in general not refer to a coordinate basis. 



orbits S in the spacetime M; notice that this excludes Kantow ski- Sachs models [19, 33]. 
We assume that these orbits S are Cauchy surfaces. Hence, M is fohated by spacehke 
Cauchy surfaces homeomorphic to S with the above symmetry property, and we let t be 
the corresponding time function. 

Up to questions about topology, Lie groups are uniquely determined by their Lie 
algebras. In order to distinguish all possible Bianchi symmetries, we therefore need to 
classify all 3-dimensional real Lie algebra; this is the Bianchi classification [33]. Let 
{^i)^2)?3} be a basis of an arbitrary 3-diixLensional real Lie algebra. The structure 
constants C°'bc with respect to this basis are defined by 

3 

c=l 

where [•, •] is the Lie bracket of the Lie algebra. In general, there exists a symmetric 
matrix {n ) and a vector with components (ai,) so that 

3 

C\c = Yl ^bcen^'' + '^bSc'' - acSb'', (4) 

e=l 

where eabc is the totally antisymmetric Levi-Civita symbol with ei23 = 1. In all of what 
follows, we restrict to the Bianchi A case, i.e., the class of uni-modular 3-dimensional 
real Lie algebras obtained by the restriction a^ = 0. Since (n ) is a symmetric matrix, 
we can assume without loss of generality a basis for which this matrix is diagonal with 
real eigenvalues ni, n2 and n^. The remaining freedom to choose the basis allows us to 
scale the numbers ni, n2 and ns, and to change the overall sign, but not to change the 
relative signs of ni, 712 and n^. In particular, if an eigenvalue is for a Lie algebra with 
respect to one basis, then there is a zero eigenvalue for all bases. The Bianchi A class 
of real 3-dimensional uni-modular Lie algebras therefore consists of 6 distinct classes, 
which are labeled as in Table 1. 

In order to apply the Bianchi classification to the construction of Bianchi A space- 
times, we introduce an orthonormal frame {60,61,62,63}. The vector field 60 is chosen 
as the future directed unit normal of the symmetry hypersurfaces Tn given by t = const 
where t is the time function above. The spatial frame vectors {61,62,63} are therefore 
tangential to each S^. As discussed in [33], it follows that 60 is a geodesic twist-free 
vector field. We can pick coordinates {t,x°') on spacetime so that {x°') are coordinates 
on each leaf Sj and cq = dt- On each T,t, we attempt to choose the spatial frame vectors 
such that [ca, Cb] = for all a, 6 = 1, 2, 3, where Ci, C2, Cs is a basis of the Lie algebra of 
Killing vector fields associated with the Bianchi symmetry. In this case, the orthonor- 
mal frame is called symmetry group invariant. Since we solve Einstein's equations as 
an initial value problem below, we can, a priori, only guarantee that the orthonormal 
frame is symmetry group invariant on the initial hypersurface Sq. It is then a matter 
of transporting the frame appropriately during the evolution in order to guarantee that 
symmetry group invariance is preserved, see below. If the orthonormal frame is symme- 
try group invariant on Sj, then {61,62,63} span a 3-dimensional real Lie algebra there 



Bianchi I 


ni = n2 = ns = 


Bianchi II 


n2 = n^ = 0, ni > 


Bianchi VIq 


ni = 0, 77,2 > 0, ns < 


Bianchi VIIq 


771 = 0, 772 > 0, 773 > 


Bianchi VIII 


ni < 0, 712 > 0, 713 > 


Bianchi IX 


77i > 0, 772 > 0, 773 > 



Table 1: The Bianchi A classes. The quantities ni, 77,2 andn^ are the eigenvalues of the matrix 
(n"'') of the tensor defined in Eq. (4). 

which is isomorphic to Lie algebra spanned by {Cij C2) Cs}- We assume that the Bianchi 
type of this algebra is one of the types in Table 1. The matrix (77"*), which we use 
in the following, is the one associated with this Lie algebra. One can show that under 
suitable assumptions below, the Einstein-scalar field equations imply that the Bianchi 
type is preserved during the evolution, i.e., if we prescribe initial data of a particular 
Bianchi type, then the corresponding solution has the same Bianchi type at all times of 
the evolution. It makes therefore sense to speak of the Bianchi type of the solution. 

From now on, tensor indices, a,b, . . . = 1, 2, 3, represent components of tensor fields 
with respect to the spatial basis vectors {ei, 62, 63}. 

2.3 The full system of equations 

Let us start with a brief summary of the derivation of the full set of equations for the 
minimally coupled self-gravitating Bianchi A scalar field case now. More details can be 
found in [33] for the very similar case of a perfect fiuid coupled to Einstein's equations 
in the Bianchi A case. The explicit form of the equations is given afterwards. With the 
choice of the frame {eo,ea} as above, we introduce the Hubble scalar H and the shear 
tensor Gab from the identity 

Ve„(eo)fe = Hhab + <yab-, 

making use of the fact that eo is a geodesic twist-free timelike vector field. Here hab is 
the induced spatial metric on the t = coTisi-surfaces. The shear tensor aab is symmetric, 
tracefree and spatial, i.e., c^i/eQ = 0. The dimension of the Hubble scalar H is (time)~^. 
We can now express all Lie brackets [cj, Cj] for 7, j = 0, . . . , 3 in terms of the objects -ff, 
0"a6; n""^ and a so far unspecified antisymmetric spatial matrix Q.ab '■= gC^eo^a, ej,) which 
determines the type of the transport of the spatial frame in time along the congruence 
tangent to cq. The Jacobi identities then imply evolution equations for each component 
of 77"*. By an "evolution equation" we mean in the following an equation which involves 
derivative terms along the timelike field eo = dt. Due to spatial homogeneity we are 
allowed to cancel all terms involving derivatives along the spatial vector fields e^ as long 
as the spatial frame vector fields are symmetry group invariant (see below). Effectively, 
all evolution equations therefore become ordinary differential equations. 

Having expressed the Lie brackets of the frame vectors as above, we can now construct 
the Riemann tensor algebraically in terms of these Lie brackets and frame derivatives 
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thereof. In particular this yields the Einstein tensor. At this stage one checks readily 
that the thus obtained expression for the Einstein tensor is a (0, 2)-tensor field which is 
symmetric if and only if the Lie brackets satisfy the Jacobi identities (as mention above) . 
In order to evaluate Einstein's equations Eqs. (1) and (2) we must express the energy 
momentum tensor appropriately (more details below). Of particular importance here is 
to assume that the scalar field respects the Bianchi symmetry. Then, the 00-component 
of these equations with respect to our frame yields the Hamiltonian constraint, the Oa- 
components imply that the commutator of the shear matrix (aab) and the matrix (n"^) 
vanishes, and the a5-components of Einstein's equations yield equations which involve 
CQ-derivatives of aab and of H. The trace of these a6-components is an evolution equation 
for H, the Raychaudhuri equation, if and only if aab is identically tracefree during the 
evolution (as we assume), and the tracefree part of the a6-components yields evolution 
equations for each component of aab- 

Notice that there is no evolution equation for the matrix 0,ab- It is therefore a 
free function which can be used to fix the remaining gauge freedom. Now, the before- 
mentioned condition that the shear matrix (aab) commutes with the matrix (n ) implies 
that initial data for which these two matrices are diagonal are no loss of generality. 
We hence assume such data in all of what follows; notice that this is consistent with 
the symmetry group invariance of the frame at the initial time. Given such data, we 
check the evolution equations above and find that the off-diagonal components of these 
matrices stay identically zero during the evolution if and only if we set Qab = (i.e., 
the orthonormal frame is parallel transported along eo). Moreover, this choice of gauge 
implies that the spatial frame is symmetry group invariant at all times of the evolution 
(and hence not only at the initial time). With this, the remaining full set of unknowns 
is H, o"i, (72, (73 (the diagonal elements of aab), ni, n2, n^ (the diagonal elements of 
Uab), the scalar field (p and the components of the orthonormal frame with respect to 
the coordinates. Since the shear tensor is tracefree, only two of the three unknowns o"i, 
o"2, C73 are relevant; we define as in [33], 

cr+ := -(a2+cr3), CT_ := — =((72-0-3). 

Now, in order to write the equations in the final form, we introduce Hubble normalized 
quantities [34] 

S±.= ^, Na.= -, 

and, for the scalar field, 

where the latter two can be interpreted as (the square roots of) the Hubble-normalized 
kinetic and potential energies, respectively, of the scalar field. These quantities are well- 
defined as long as H does never become zero during the evolution. We assume here 
that H > initially and hence that the universe is expanding initially. Below we show 
that for all Bianchi A models, possibly except for Bianchi IX, it follows that H > 



during the whole evolution and that hence, possibly except for Bianchi IX, these Hubble 
normalized variables are well-defined. 

Now, dividing the Hamiltonian constraint mentioned above by H^ yields the Fried- 
mann equation 

l = ^l + ^l+x^ + y^ + K. (6) 

Here, 

^■=-■^ = -[^(^1+ ^2 + Ni - 2{NiN2 + N2NS + N^Ni)) , (7) 

where ^i? is the spatial Ricci scalar. The Raychaudhuri equation mentioned above 
becomes 

eo{H) = -{l+q)H^, (8) 

with the deceleration scalar 

q := 2S2 + 2x^ - y^. (9) 

Here and in the following, we use the shorthand notation S^ := S^ + S^. Notice that 

the deceleration scalar q measures the acceleration of the expansion of the cosmological 

model. Indeed, if g > 0, then the expansion is decelerated (hence becomes slower), while 

if g < 0, the expansion is accelerated (hence becomes faster). In particular, a change of 

sign from negative to positive during the evolution signals a graceful exit. 

In order to write down the other equations above, we introduce a new time coordinate 

r, the so-called Hubble time, by 

dT _ 

Tt-^' 

and we denote derivatives with respect to r by a prime '. Then, by dividing the before- 
mentioned evolution equations for the shear variables by H^, we find 



where 

N^)) , (10) 

(11) 

Notice that 5"+ and S- are linear combinations of the eigenvalues of the Hubble nor- 
malized tracefree part of the spatial Ricci tensor. In the same way we proceed with the 
before- mentioned evolution equations for ni, 722 and n^ and obtain 





E'± = -(2-(7)S±-5±, 


'+ 


= ^{(N2-N^f-Ni{2Ni-N2 




= -J-(N^ - N2){Ni - N2 - N^) 



iV^ = (Q + 2S+ + 2V3S_)iV2, 
N^ = {q + 2S+ - 2\/3S_)iV3. 



Next we extract equations from Eq. (3). Under the assumption of Bianchi symmetry, 
it reads 

eo(eo(</')) + 3/feo(0) + y'(</.) = O. 

This equation, together with the definition of x and y above, gives 

x' =x{q-2)-F{<i))y\ 
y' = F{(j))xy + y{l + q), 



where we define 



-<^'-\/|^- 



The evolution equation for cp is therefore simply 

(p' = vGx. 

Since we can rewrite Eq. (8) as 

H' = -il + q)H, (13) 

and hence interpret this as an evolution equation for H and since we can derive evolu- 
tion equations for the coordinate components of the orthonormal frame easily from the 
condition that the connection is torsion-free (we shall refrain from writing these equa- 
tions down now), we have therefore now succeeded in obtaining the complete system of 
equations for all unknown variables. In all of what follows, we shall focus on the main 
core of the evolution equations 

S'± = -(2-(?)E±-5±, (14) 

iV( = (g-4S+)iVi, (15) 

iV^ = (g + 2S+ + 2^/3S_)iV2, (16) 

iV^ = ((7 + 2S+-2^/3S_)iV3, (17) 

x'=x{q-2)-F{^)y^ (18) 

y' = F{<P)xy + y{l + q), (19) 

0' = VQx, (20) 

where, 5*+, S- and q are given by (10), (11) and (9) and where the function F((j)) is 
considered as a given function of (p defined in Eq. (12). This set of evolution equations is 
a dynamical system whose state space is spanned by (S+, S_, Ni, N2, N3, x, y, (p) G M^. 
It is subject to the Friedmann constraint Eq. (6). The quantity 

G:=l-^l-T.l-x^-y^ -K, (21) 

is a measure of the violation of this constraint. Suppose we prescribe data for the initial 
value problem of the evolution equations above which possibly violate the constraint. 



i.e., G does not necessarily vanish initially, and then consider the corresponding solution 
of the evolution equations. Then it is straightforward to derive an evolution equation 
for G, namely 

G' = 2qG, (22) 

the subsidiary system. This implies that if we prescribe initial data which satisfy the 
constraint, i.e., G = initially, and then determine the corresponding solution of the 
evolution equations, then the constraint is satisfied identically for all times. 

After the main set of evolution equations Eqs. (14) - (20) above has been solved, 
Eq. (13) can be used to determine H as a function of r. Moreover, the coordinate 
components of the orthonormal frame can be computed as a function of r as described 
above. 

In our application, the scalar field (p often approaches infinity rapidly. It is then useful 
to replace the variable 4> hy ip := !/(/>. With this we obtain the following alternative 
evolution equations 

x' = x(q-2)-F{i;)y^, (23) 

y' = F{i,)xy + y{l + q), (24) 

ip' = -yflx'4?. (25) 

Here, we use a slightly sloppy notation where the function F, which was defined as a 
function of </> above, is now considered as a function of ^ = l/i?i>, i.e., 

3 Numerical implementation 

3.1 Initial data 

For solving the equations numerically, we are interested in initial data that satisfy two 
conditions: (i) They must satisfy the constraint Eq. (6) and, (ii) they should represent 
an inflationary epoch, i.e., the quantity q in Eq. (9) should be negative initially. Finally, 
(iii), the Bianchi type of the model is also determined from the initial data and hence 
the data should be chosen in a particular way. 

We proceed as follows. Eqs. (6) and (9) can be solved for x^ and y^ and we get 

x2 = (1-3^2 -K + g)/3, y^ = (2-2K-q)l'i. 

This implies the restrictions 

1 - 3^2 - K + g > 0, 2 - 2K - g > 0, 

since we shall always demand that y > (see the discussion in Section 4.2) and hence 

-iY? + K-\<q< 2(\-K). 
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This inequality can be satisfied for q < if and only if 

3^^ + K < 1. 

holds. Notice here that 2(1 — K) > as a consequence of the Friedmann constraint. 
Based on these insights, we adopt the following strategy to choose inflationary initial 
data for our numerical evolutions: 

1. Choose arbitrary numbers for the initial data of A^^i, A'^2 and N3 (compatible with 
a specified Bianchi class) and for S+ and S_ so that 

3S2 + K< 1. 

2. Choose an arbitrary number q in the interval 

3i;2 + /C - 1 < g < 0. 

3. Choose the initial data for x and y as 

X = ±V(l-3S2-i^ + g)/3, y = V(2 -2K- q)/3. 

4. Choose an arbitrary positive number for the initial data of the scalar field^ (j). 
In many situations, we shall pick the free parts of the initial data randomly. 

3.2 Numerical implementation and constraint damping 

In order to solve the ordinary differential evolution equations numerically, we have used 
the 4th-order Runge Kutta method (non-adaptive) and the Dormand-Prince method 
(adaptive time stepping). 

We made several numerical experiments with this setup and found quickly that the 
evolution equations in the form above have a serious numerical problem. In Figs. 1 and 2 
we see that, while the constraint violation quantity G (defined in Eq. (21)) is numerically 
well-behaved as long as q < 0, it grows rapidly during epochs with q > 0, i.e., when 
the expansion is decelerated. Eventually the numerical solutions are therefore rendered 
useless after short evolution times when G becomes of order unity. This observation 
is in agreement with the subsidiary system Eq. (22) of the evolution equations. If G 
is not exactly zero initially, we can conclude that it decays exponentially while q < 
(when the expansion is accelerated) and increases exponentially while q > (during 
decelerated epochs). At the beginning of the numerical calculations, q is negative and 
hence, the constraint violations, which are non-zero initially, do not grow; the fact that 
they do not decay, as it is suggested by Eq. (22), is caused by numerical discretization 
errors since the evolution equations are solved numerically and hence are satisfied only 
approximately. 



In Section 4.5 we discuss further restrictions for the choice of ( 
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Fig. 1: Deceleration scalar q (BI). 
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Fig. 2: Constraint violations G using the orig- 
inal evolution equations (BI). 



Before we present our solution to this problem, let us explain why the violation of 
the constraint at the initial time in Fig. 2 is '^ 10^^^ , and not, as it may be expected, 
at the level of machine precision ~ 10~^^. The reason is simply that in all of our figures, 
we do not plot the initial data themselves (which correspond to r = 0), but rather the 
solution from the first numerical evolution step onwards. 

In order to solve the instability of the constraint hypersurface of our dynamical 
system and hence to be able to make reliable long-term numerical computations, we 
propose to solve a system of modified evolution equations which agrees with the original 
one precisely if the constraint violations vanish and which is obtained by adding suitable 
constraint damping terms to the evolution equations. These modified evolution equations 
are 



s'+ = 


= -{2-q)i:±-S+ + hG, 


s'_ = 


= -{2-q)^±-S-+f2G, 


N[ - 


= (g-4S+)iVi + /3, 


N', - 


= (g + 2S+ + 2\/3S_)A^2 + /4G' 


N', - 


= (g + 2S+-2\/3S_)Af3 + /5G 


x' -- 


= a,(g _ 2) - FCV^y + /eG, 


y' - 


= F{^p)xy + y{l + q) + frG, 


V.' = 


= —\/6xip . 



Comparing this to the original "unmodified" evolution equations Eqs. (14) - (20) (or 
more precisely (23) - (25)), it becomes obvious that this modified system agrees with 
the original one precisely if G = 0. Hence any solution of the original system, which 
satisfies the constraints, is also a solution of this modified system and vice versa. The 
functions fi, . . . , fi are so far unspecified. The subsidiary system for these modified 
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Fig. 3: Constraint violations G using the modified evolution equations (BI). 



evolution equations is 

G' = -2S+S+ - 2S_S'_ - 2xx' - 2yy - K', 
which, using the evolution equations above, motivates us to choose 



/5 = -(A^i + iV2)/2, /6 = 2X, 



/t = (k - l)y. 



-(iVi + iV3)/2, 



for some constant k > 0. Certainly, this choice for /i,. . . ,/7 is not unique, but it has the 
nice consequence that the subsidiary system becomes 



G' = -((iVf + Ar| + nI)/6 + 2Ky2)G. 



(27) 



In particular, the factor in front of G on the right hand side is negative semi-definite and 
hence, the constraint violations are expected to be numerically stable irrespective of the 
sign oi q. In all of what follows we choose k = \. This theoretical prediction is confirmed 
numerically in Fig. 3 which shows numerical evolutions for the modified system. We see 
that the constraint violations are now stable during the whole evolution and are driven 
towards numerical double precision round off errors (~ 10^^^). We can now hope that 
we are able to do reliable long-term numerical simulations. 

The above examples are Bianchi I solutions. We find that this way of damping 
constraint violations works well in all Bianchi cases as long as the numerical resolution 
is sufficient. Indeed, since some of the unknowns are unbounded in several Bianchi 
cases (more details on this later), keeping a sufficient numerical resolution can become 
a challenge in practice (cf. e.g. Fig. 24). 

A possible alternative to introducing constraint damping terms to the "free" evo- 
lution equations is to implement a "constrained evolution" scheme. Here, instead of 
determining the solution for all variables from the evolution equations, we would single 
out one variable and instead determine it algebraically from the Friedmann constraint. 
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By this, we would guarantee that the constraint violations are zero (numerically). In 
order to demonstrate, however, that the resulting numerical solution satisfies all of Ein- 
stein's equations, we would then need to check if the one evolution equation, which has 
been neglected in the first step, is satisfied. We hence arrive at the same problem as 
in our implementation, namely, that there is a-priori no guaranty that this additional 
equation is satisfied well by the numerical solution. The origin of this problem is that 
Einstein's equations are overdetermined and this is the same for "free evolution" (as in 
our case) and such "constrained evolution" schemes. 

4 Analysis of the graceful exit problem 

4.1 Basic strategy 

We restrict to expanding cosmological models characterized hy H > where H is the 
Hubble scalar. The acceleration of the cosmological model is described by the decel- 
eration parameter q in Eq. (9) which is, by convention, positive when the expansion 
is decelerated and negative when it is accelerated. In the following we shall prescribe 
initial data for the dynamical system above which represent the early universe in an 
inflationary regime, i.e., for which q < initially. If the corresponding solution has the 
property that q < for all times, then the inflation era is eternal. If, however, q changes 
its sign after a finite time and stays positive, then there is a graceful exit from inflation. 
The main idea which we follow in this paper is to identify certain future attractor so- 
lutions in the decelerated regime - in many cases these are stationary solutions of our 
dynamical system. Notice that we use the term "attractor" sometimes in a slightly loose 
sense in this paper. In the strict mathematical sense, an attractor is, in particular, a 
subset of the state space which is approached by generic orbits (the precise definition 
can be found in [33]). In the Bianchi cases VIIq and VIII below, however, we refer to 
a particular generic asymptotic behavior as an attractor even though this asymptotic 
behavior cannot be associated directly with a subset of the state space. 

Stationary solutions for homogeneous and isotropic scalar field models were studied 
before, see for instance [11, 20]. If future attractors in the decelerated regime exist then 
our solutions must have the property that the initial inflationary epoch is finite and is 
succeeded by a decelerated regime after a graceful exit. 

4.2 Basic consequences for the evolution 

Here we summarize some basic results about the dynamics of our models. Let us assume 
from now on that the potential y is a strictly positive function of (p. 
Consider the Friedmann equation Eq. (6), which can be rewritten as 

H^ = al + a^_, + (eo((/>)) Ve + ^(0)/3 - ^R/6, 

after multiplication with H^. Let us assume that H > at the initial time. Hence H 
can become non-positive during the evolution only if the right side of this equation is 
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allowed to become zero. Under our assumptions, however, this is only possible if ^R can 
be positive. It is a standard result, see for instance [33], which can be checked directly 
for instance from Eq. (7), that ^i? is non-positive for all Bianchi A models possibly 
except for Bianchi IX. For all other Bianchi A cases, we have in fact that < K < 1. By 
excluding Bianchi IX for the rest of this paper, we therefore make sure that H is positive 
during the whole evolution and that hence our cosmological models do not recollapse. In 
particular, this means that the Hubble normalized variables, which we have introduced 
in the previous section, are well-defined during the whole evolution. 

We can also show that our cosmological models must isotropize and spatial curvature 
must decay during inflation (i.e., while q < 0). To this end, define a new quantity 

S — l-x^-y"^. (28) 

The Friedmann constraint implies that S = Ti'^ + K and < S" < 1 (if we exclude the 
Bianchi IX case). The evolution equations (14) - (20) yield an evolution equation for S 

S' = 2qS - 4S2 < 2qS. (29) 

Thus, 5, and therefore E^ and K, decrease rapidly during inflation {q < 0). With some- 
what more work not discussed here, one can show that in fact the full spatial curvature 
tensor decays during inflation. If inflation lasts forever, then this is a realization of the 
cosmic no-hair conjecture. We are here, however, rather interested in situations when 
inflation does not last forever and there is a graceful exit. If inflation is sufficiently long, 
then our argument here suggests that anisotropies and the size of the spatial curvature 
are extremely small by the time of the graceful exit. One of the questions which is of 
interest for us here is then: What happens when inflation is over? Do anisotropies and 
spatial curvature stay small or do they grow again? 

Another important immediate consequence from the fact that K > for all Bianchi A 
cases (possibly except for Bianchi IX) is that the constraint Eq. (6) implies the bounded- 
ness of the variables Ti± , x and y. From the definition of K in Eq. (7) and the conventions 
for the Bianchi cases in Table 1, we find that the remaining quantities have the follow- 
ing properties. For Bianchi I, the variables Na are identically zero and hence trivially 
bounded. For Bianchi II, the non-zero quantity A'^i must be bounded. For Bianchi VIq, 
the two non-zero quantities A'^2 and A^3 (with opposite signs) are bounded. For Bianchi 
VIIq, however, the two non-zero quantities A'^2 and N^ (with equal signs) are not neces- 
sarily bounded. For Bianchi VIII, it follows that A^i must be bounded, but A'^2 and A'^2 
can be unbounded. 

We also mention here the following basic inequalities and extreme cases for the 
deceleration scalar q. The formula Eq. (9) for q can be rewritten as 

q = 2 - 3y2 - 2K, (30) 

using the Friedmann constraint. Hence we have q < 2 for all Bianchi A cases (possibly 
except for Bianchi IX) and g = 2 if and only if i^ = and y = 0. We can therefore 
typically assume that q < 2. Since y^ + K < 1 (from the Friedmann constraint), it also 
follows that q > —1- The case q = —1 corresponds to de-Sitter like exponential inflation. 

15 



4.3 Known results and our choice of the scalar field potential 

Let us now list the main relevant rigorous results for minimally coupled Bianchi A scalar 
field models for various classes of potentials. We shall use those insights in order to 
choose a concrete class of potentials for our numerical investigations for which solutions 
with graceful exits from inflation are possible. 

The first relevant result has been obtained by Rendall [26]. If the potential has a 
positive lower bound (and certain further technical conditions are satisfied), then the 
solutions undergo eternal exponential de-Sitter like expansion asymptotically. Hence 
potentials with a positive lower bound do not allow solutions with a graceful exit from 
inflation. Since we want to avoid the complication of chaotic inflation associated with 
potentials with a zero local minimum [23, 2], we must therefore focus on potentials which 
approach zero in the limit (p — t- oo. This could be an exponential potential [15, 21, 3], 
in which case it was shown that eternal power-law expansion occurs if the exponent is 
not too negative. We shall return to this potential later (but it will not be our main 
choice). Another possibility is an "intermediate" potential. The investigations in [27] 
(see in particular Theorem 4 there) yield that one gets eternal "intermediate inflation" 
if, in particular, lim0_j.oo V^'((/')/V^('/') -^ 0. See also [5]. We must therefore consider 
potentials with the property liuifp^oo V{4')/^{4') < in order to make solutions with 
graceful exits from inflation possible. 

Now the concrete model, which we have decided to study in this paper, is the potential 
introduced by Albrecht-Skordis [1] in the context of FRW cosmological models with dark 
energy in string theory: 

y(0) = e-^i^C2 + </.2), (31) 

where ci > 0, C2 > are constants. In particular, the function -F(V') (see Eq. (26) where 
(p has been replaced hy (p = l/ip) has the form 

Its limit 

hm F(V') = - J|ci =: Fo < 0, (33) 

i/)-5>0 V 2 

and hence one of the sufficient conditions for eternal intermediate infiation in [27] men- 
tioned above is in general violated. 

As shown in Fig. 4, the potential Eq. (31) gives rise to three different cases. According 
to the idea that the scalar field essentially "roles down the potential hill" during the 
evolution, we expect that the scalar field either goes to infinity asymptotically (in which 
case the discussion above suggests that a graceful exit from inflation is possible), or is 
stuck in the local minimum of the potential where we then expect eternal de-Sitter like 
inflation according to [26]. Since we are interested in the graceful exit problem, we shall 
start our investigations for cases when ci and C2 are such that the potential is strictly 
monotonic. One shows easily that this is the case when cfc2 > 1. Later in Section 4.5, 
we present numerical investigations for other choices of the parameters of the potential, 
in particular when the potential has a "valley" . 
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Fig. 4: The three main cases for the scalar field potential V {({)). 

4.4 The graceful exit problem for the monotonic case of the potential 

4.4.1 Bianchi I (iVi = A^2 = A^3 = 0) 

Let us start our discussion with the simplest case: Bianchi I. This is characterized by 
the condition that A^i, A^2 and A^3 vanish identically which implies spatial flatness. 

Future attractors. We first identify a Bianchi I stationary solution of the dynamical 
system, i.e., a solution of the evolution equations Eqs. (14) - (17), (23) - (25), and of 
the constraint Eq. (6) with the property 



S'± = N[ = N!, 



AT^ 



y 



^' = 0. 



We determine this stationary solution under the condition ip = (which represent the 
hmit </>—)• oo when the scalar field has "completely rolled down the potential hill" 
asymptotically) and for y > 0. In the Bianchi I case, the condition S^ = yields 
S-t = (see Eq. (14)) since 5-i- = and since g < 2 for y > (see Section 4.2). Notice 
that y > is suggested by the numerical solutions below (as opposed to y = 0). We 
can use the Friedmann constraint to express y in terms of x, i.e., y = \/l — x^ . The 
evolution equation for x with x' = can now be solved for x. This algebraic equation 
yields two solutions of which only one is relevant. The evolution equation for y with 
and y = \/l — x^ is then satisfied identically. We get 



y 



E± = 0, iVi = N2 
with Fo given by Eq. (33), i.e., 



N:^ 



0, X 



9 ' ^ 



(34) 



Fn 



Cl. 



Hence, for every fixed ci < \/6j we have a stationary solution of type Bianchi I. Notice 
here that this is strictly speaking not an actual solution of the field equations because 
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^ = (corresponding to i;^ = oo) lies outside the allowed range. This stationary solution, 
however, can be interpreted meaningfully as the limit of actual solutions when — )• cxd. 
We shall always understand our stationary solutions in this way. Eq. (9) yields the 
corresponding value of q: 

'-1 



c?-2 



2 
It follows that q > Q (i.e., the stationary solution is in the decelerated regime) if and 
only if 

V2 < ci < Ve. (35) 

We claim now that generic Bianchi I solutions with our choice of the potential and for 
ci compatible with Eq. (35) approach this stationary solution in the future asymptoti- 
cally and hence that it is a future attractor. If this is the case the solution is decelerated 
and becomes isotropic and spatially flat. In particular, for evolutions starting from infla- 
tionary initial data, there is a graceful exit after a finite evolution time and the process 
of isotropization and approach towards spatial flatness, which occurs during inflation 
according to our general discussion in Section 4.2, continues after inflation. 

If this stationary solution is really a future attractor for Bianchi I models, then it is a 
necessary condition that it is future non-linear stable. According to the standard theory 
of dynamical systems, we linearize the evolution equations for S-t, x, y and ip with 
A'^i, A^2) -^3 identically zero around the stationary point. We find that the corresponding 
evolution matrix has the eigenvalues (cf — 6)/2 (three times repeated), c\ — 2 (single) and 
(single). The stationary point therefore has a 3-dimensional future stable subspace, 
a 1-dimensional future unstable subspace and a 1-dimensional center subspace under 
the restriction Eq. (35). We check that the 1-dimensional future unstable subspace is 
transverse to the constraint hypersurface and hence does not occur for solutions which 
satisfy the constraint. 

In order to understand the meaning of the center subspace, we apply the center 
manifold theory; see for example the summary and references in [29]. For this discussion, 
we throw out the evolution equation for y and replace all occurrences of y in the other 
evolution equations by the expression implied when the Friedmann constraint is solved 
for y (for y > 0). The center subspace is then the 1-dimensional subspace of the tangent 
space at the stationary point given by S-t = 0, x — ci/\/6 + y^2/3'(/' = 0. This means 
that ij) \s a. regular local coordinate of the corresponding center manifold (which exists 
but is not necessarily unique according to the center manifold theory). The evolution 
equations yield that the center manifold can be approximated in terms of this local 
coordinate as 




S± = O('0), 2;-ci/V6+ V2/3^ = -W- ^^ +0(V'), at ^ = 



Also higher orders in -0 can be computed easily. We can therefore replace x in the 
evolution equation (25) for ^ by this expansion and find a new evolution equation for ^ 
which only involves V'^ 

^' = -ciV'^ + 0(V;3^. 
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Fig. 5: Scalar field kinetic energy x (BI). 
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Fig. 8: Scalar field ip = 1/0 (BI). 



It then follows that all solutions iP{t) of this differential equation approach the value 
of the stationary point, at least if the initial value of -0 is sufficiently small. Hence, the 
center manifold is future stable. We compute that the solutions on the center manifold 
have the asymptotic expansions 



V'(r) = -T i + 0,r-2 + ..., x(^t) = ^ 



ci 



2 1 

3 ci 



T-' + .. 



at r — 7- oo, where 0* is a free datum which also appears in higher terms in the expansion 
of x(r). 

Numerical studies. We have therefore demonstrated that the stationary solution is 
future stable in the Bianchi I class. This does of course not mean that it is a future 
attractor. We shall now check the attractor property numerically, but provide no proof. 
For the numerical computations, we choose ci = 1.430, C2 = 1.0, which is in agreement 
with Eq. (35). Then, the stationary solution is T,^ = 0, S_ = 0, x = 0.583, y = 0.811 and 
q = 0.022. In Figs. 5 and 6, we demonstrate that the solutions corresponding to different 
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inflationary data indeed approach these values asymptotically in the future. We point 
out that it is crucial here to use the modified evolution equations discussed in Section 3.2 
in order to get reliable numerical results; in particular close to the stationary point where 
there exists an unstable constraint violating mode. In any case, we have hence indeed 
constructed models for which inflation lasts only for a flnite time and a graceful exit 
into a permanent decelerated epoch follows. Fig. 7 shows that the numerical result is 
consistent with the exponential decay of the shear variables with exponent (cf — 6)/2 
(the bold straight line in Fig. 7) close to the stationary point as implied by the stability 
analysis above. Moreover, the 1/t-decay of the function ifj in Fig. 8 as predicted from 
the center manifold analysis close to the stationary point is consistent. 

As a further representation of the dynamics in the Bianchi I case, consider Fig. 9. 
In this flgure, we plot the projection of the evolution orbits for three different initial 
data to the x-y-plane in the state space for the same case ci = 1.43, ci = 1.0 as before. 
According to the Friedmann constraint, the dynamics must take place in the interior of 
the circle x^ + y^ = 1 (denoted by the dotted ellipsoidal curve in the figure), and on 
this circle if and only if Tj± = (i.e., the isotropic case). According to Eq. (30), the 
horizontal line y = y^2/3 ~ 0.82 in the figure corresponds to g = 0, while everything 
above this line yields g < (inflation) and everything below is g > (deceleration). As 
we see, the three choices of initial data start in the inflationary regime (i.e., above the 
q = 0-line). In accordance with the discussion of the quantity S defined in Eq. (28) 
in Section 4.2, the orbits approach the circle x^ + y"^ = 1 rapidly and hence become 
isotropic very quickly during infiation. In fact, all orbits become close to the unit circle 
so quickly that some of their parts can hardly be distinguished in the figure. The orbits 
follow (as shown in the figure) the unit circle in such a way that eventually they all cross 
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the q = 0-line (graceful exit) and then approach the stationary point discussed before. 

Results. All this suggests that the stationary solution above is indeed a future at- 
tractor for Bianchi I models in the decelerated regime if v2 < ci < \/6. Hence if we 
start with initial data in the inflationary regime then a graceful exit must occur under 
generic conditions. Since the stationary solution is isotropic, anisotropies continue to 
decay even after inflation. We stress that only the limit lim^_!.oo ^'(0)/^('/') is relevant 
for the existence and main stability properties of the stationary point above. Hence it is 
conceivable that a similar dynamics can be found for more general classes of potentials. 
A particularly important example is the exponential potential studied in [21]. For this 
potential only the 0^-term is missing in Eq. (31). In particular, the parameter ci has the 
same meaning, and one finds that the limit in Eq. (33) is the same. We therefore con- 
jecture that generic inflationary Bianchi I with an exponential potential have a graceful 
exit after a finite time under the same conditions as here, namely if y/2 < c\ < -\/6- 
This would be remarkably consistent with the main result of [21], namely, that there 
is eternal power-law inflation for < ci < v2. As a final remark, let us comment on 
the asymptotic behavior of the Hubble scalar H for our models. Since q approaches a 
stationary positive value asymptotically, H decays towards zero exponentially according 
to Eq. (13). 

4.4.2 Bianchi II {Ni > 0, N2 = N3 = 0) 

Future attractors. In the Bianchi II case, we require the variable A'^i to be positive 
while A'^2 and A''3 vanish. It is conceivable that the Bianchi I stationary solution in the 
previous section plays the role of a future asymptotic end state. However, when we 
repeat the stability analysis, but this time allow for perturbations with A^i 7^ 0, it turns 
out that the subspace of the tangent space generated by A^i is unstable. Hence, the 
above Bianchi I stationary point is not relevant in the Bianchi II class. 

We therefore construct now a genuinely Bianchi II stationary solution with y > 
and ip = as before. As in the Bianchi I case, the condition S'_ = implies that S_ = 
since S- = 0, see Eq. (11). But, since 5+ does not necessary vanish, see Eq. (10), the 
condition S^ = yields that E_|_ can be non-zero. The resulting algebraic system of 
equations enables us to find the following stationary solution of the Bianchi II dynamical 
system 



S_=0, A^i= ^' ,, "i; ° -, N2 = N3 = 
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As before, Fq can be expressed in terms of ci using Eq. (33). The solution is real, and 
satisfies y > and g > if and only if 



\/2 < ci < \/8. 



(36) 



The same analysis as in the Bianchi I case reveals that this stationary point is future 
stable. 

Numerical studies. We claim that this stationary solution is a future attractor for 
Bianchi II solutions, and support this claim by numerical studies. For this case, we 
choose ci = 1.450, C2 = 1.0 which meets Eq. (36). The corresponding stationary values 
are S+ = 0.011, S_ = 0, A^i = 0.257, x = 0.588, y = 0.804, and q = 0.045. In Figs. 10, 
11 and 12, we demonstrate for various inflationary initial data that the evolution indeed 
approaches these values asymptotically. In particular, a graceful exit from inflation 
always happens. 
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Results. We have therefore identified a future attractor of type Bianchi II in the 
decelerated regime if \/2 < ci < v8. Hence under this condition, generic initially 
inflationary solutions have a graceful exit from inflation very similar to the Bianchi 
I case. The main difference to the Bianchi I case is the asymptotic behavior of the 
spatial curvature and anisotropy variables. Because the values for S_|_ and A^i of the 
stationary attractor are non-zero, isotropization and decay of spatial curvature stops at 
the time of the graceful exit. Even more so, as we see in Figs. 10 and 11, these quantities 
can grow again rapidly after inflation. In particular, this brings doubt to some basic 
assumptions of the standard model. Namely, the apparent isotropy of our universe may 
in general not be satisfactorily explained by a finite phase of rapid expansion. As our 
example shows it may in fact be possible that, without an additional so far unknown 
mechanism, anisotropics and spatial curvature, which are presumably very small by the 
end of inflation, grow again when inflation is over. 

4.4.3 Bianchi VIq (iVi = 0, A^2 > 0, N3 < 0) 

Future attractors. For the Bianchi VIq case, we assume A'^i = 0, and that A'^2 and 
N3 have opposite signs. Again the Bianchi I stationary solution above could be a future 
asymptotic end state for general Bianchi VIq orbits. In the same way as in Bianchi II, 
however, it is future unstable in Bianchi VIq. The Bianchi II stationary solution above 
could also be a future asymptotic end state after applying a rotation to the orthonormal 
frame. However, also this is unstable in Bianchi VIq. 

We therefore attempt to construct a Bianchi VIq stationary solution with q > and 
y > as before. The condition N2 = N^ = 0, but N2,N3 ^ 0, implies that q = — 2S+ 
and E_ = 0, see Eq. (16) and (17). Then Eq. (14) together with Eq. (11) yields that 
A'^2 = —-^3. Using these relations, we find that following result 

6 + F^ 6 + F(f 

X = Z-TT, V = —77, "0 = 0, 



where 



,^^1-^ 



F2 + 6- 
Expressing Fq in terms of ci as before, we find that the only restriction is 



ci 



> V2. (37) 



Results. Since the results are so similar to the Bianchi II case, except that in addition 
here the modulus of A'^2 and N3 approaches the same value asymptotically while A^i = 0, 
we do not show any numerical results now. We come to the same conclusions as in the 
Bianchi II case that for generic infiationary Bianchi VIq initial data there is a graceful 
exit if ci > v2- Again the future attractor is neither isotropic nor spatially flat. 
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4.4.4 Bianchi VIIq (iVi =0, N2,N3> 0) 

We proceed with the Bianchi VIIq case where A'^i is assumed to vanish and A'^2 and N^ are 
both positive. At a first glance, the situation seems to be very similar as in the Bianchi 
VIq case. However, as N2 and N^ have the same sign now, the Friedmann constraint 
does not imply boundedness of A''2 and N3 . 

Future attractors. The numerical studies which we present in more detail below 
suggest that for generic inflationary Bianchi VIIq initial data, A''2 and N3 are unbounded. 
This means that we are now not looking for stationary solutions in order to describe the 
future asymptotics. However, all quantities S-t, x, y and ^ seem to approach stationary 
values (as it was the case for all previous Bianchi cases), namely that 

S-t — )• 0, X— )-x=K, y— )-y*, tp ^ 0, for r — )• 00, (38) 

where the values x=k and y=K are universal and only depend on ci. We claim now that for 
every choice of the parameters ci and C2 (in a certain range) , we can construct a future 
attractor for Bianchi VHq solutions in the decelerated regime as follows. According to 
Eqs. (9) and (38), the deceleration scalar q must approach a universal value q^: = 2x1— y^. 
Eqs. (16) and (17) then imply that 

N2{t) ^ N2,e'"\ Nsir) ^ Ns.e'''^ , (39) 

for large r, where A''2* and N^^, are some positive constants. On the other hand, Eqs. (14) 
with Si -> and T,'^ -^ imply that S+,S- -^ 0. According to Eqs. (10) and (11) 
these quantities read in the Bianchi VIIq case: 

S+ = l{N2-N,f, S. = ^{Ni-Ni). (40) 

This implies in particular that A^2* = ^3* (but also gives a restriction on the next order 
of the asymptotic expansion for A'^2 and N3 which we do not discuss here in detail now). 
The quantity K, which according to Eq. (7) is K = {N2 — N^)"^ /2 in the Bianchi VIIq 
case, must also approach zero. We can therefore consider Eq. (6), (18) and (19) easily 
in the limit r — )• 00, and find algebraic equations for x=k and y-t which are the same as in 
the Bianchi I case. The solution is 



^/6' 



In particular, if 



(41) 



\/2 < ci < Ve, (42) 



we claim that generic initially inflationary Bianchi VIIq solutions approach the asymp- 
totic behavior given by Eqs. (38), (39) with A^2* = -^3* and (41) in the future. 
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Fig. 13: Deceleration scalar q, (BVIIq). 
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Fig. 14: Oscillatory decay of \S+\ (BVIIq). 
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Fig. 15: Oscillatory decay o/|5_| (BVIIq). Fig. 16: Late time oscillations of S- (BVIIq). 

Numerical studies. For our numerical studies, we choose ci = 1.45 and C2 = 1.0. 
This implies that x^: ^ 0.59, y* ?» 0.81 and g* ss 0.05 according to Eq. (41). Indeed, 
the variables E-t, x, y and ip approach these stationary values generically, and hence in 
particular q (see Fig. 13). We therefore find that our initially inflationary Bianchi VIIq 
solutions have graceful exits from inflation. Of major interest for the Bianchi VIIq case is 
the behavior of the variables A'^2 and N^ which are supposed to behave like Eq. (39) with 
^2* = N^* so that 5_|_, 5_ — )■ 0. In Figs. 14 and 15 we see that this decay towards zero 
indeed happens exponentially. However, the decay is not monotonic and the frequency 
of the oscillations increases as time progresses. While 5+ is non-negative, the quantity 
S- can have any sign, cf. Fig. 16. The evolution equations for T,± therefore suggest that 
also S-t should be oscillatory. The plots in Figs. 17 and 18 confirm that while S-|_ decays 
monotonically with the same rate as in the Bianchi I case (represented by the bold line), 
the quantity S_ decays with a smaller rate than in the Bianchi I case and is oscillatory 
with increasing frequency. 



25 




-2 


^ 


-^- 


^\^^ 


-4 
-6 
-8 


"\\ ; 


10 
12 




Initial data 

id.1 

id.2 

id.3 


\ \ 


14 






\ , \ ; 






_^4— ^ 


-5 


^--ip 


\ 




10 




initial data 


\ 


m 






id.1 


\ 


11. 






id.2 


\ ■ 


n. 






id.3 


\ 


?ft 


15 






\ 






:;"|,. 



200 220 240 260 280 300 



200 220 240 260 280 300 



Fig. 17: Decay of T.+ (BVIIq). 



Fig. 18: Oscillatory decay ofT.^ (BVIIq). 



Results The Bianchi VIIq case differs significantly from the previous Bianchi cases due 
to the unboundedness of the quantities N2 and A^3 . It is particularly interesting to note 
that the situation is not the same in pure vacuum where A'^2 and N3 have been proven 
to be bounded for generic Bianchi VIIq solutions; see Theorem 1.1 in [30]. In our case 
here, A''2 and N^ grow exponentially as in Eq. (39) where q^: is given in Eq. (41) such 
that 5+ and 5_ given by Eq. (40) go to zero. This leads to a subtle oscillatory behavior 
which is passed on as oscillations for the otherwise exponentially decaying shear variables 
Ti±. This dynamics should be studied in more detail in future work, in particular, the 
dependency of the oscillation frequency on time and the decay rate of S_ in Fig. 18. 
Nevertheless, our main claim that generic initially inflationary Bianchi VIIq solutions 
have a graceful exit from inflation (if ci is chosen appropriately), has been confirmed. 
We have demonstrated that the variables T,±, x, y and ip approach stationary values 
similar to previous Bianchi cases. Since the future attractor is isotropic in the Bianchi 
VIIq case, anisotropies (represented by S-t) and the spatial curvature (represented by 
the variables 5*+, S- and K) continue to decay even after infiation. 

4.4.5 Bianchi VIII (A^i < 0, iV2,A^3 > 0) 

The most complicated case, which we study in this paper, is Bianchi VIII where we have 
the negative variable Ni in addition to the positive variables A'^2 and N^. Similar to 
Bianchi VIIq, the quantities N2 and N^ are not bounded by the Friedmann constraint, 
but S-t, X, y and A'^i are. 

Future attractors. The numerical studies presented below suggest that all quantities 
S-t, X, y, Ni and tp approach stationary values during the evolution, but that A'^2 and 
N^ are unbounded. Hence the situation appears to be quite similar to the Bianchi VIIq 
case. The numerical computations suggest that 



j+*, 



0, A^i -^ 0, 



X*, y-j-y*, tp-^0, for r 



00, 



(43) 
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where the values x*, y* and £+=„ are universal and only depend on the parameter ci of 
the potential. We notice that in particular S+ does in general not approach zero in the 
Bianchi VIII case (in contrast to Bianchi VIIq). According to Eq. (9), the deceleration 
scalar q approaches the universal value 

q, = 2S^, + 2x1 - yl (44) 

Eqs. (15) - (17) imply that 

A^i(r)^iVi.e(^*-^^+*)", A^2(t) ^ A^2*e(''*+2^+*)", N^{t) ^ N^J'i'+^^+'^\ (45) 

for large r with constants A''i* < 0, A'^2* > and A'^s* > 0. Eqs. (10) and (11) read 

S+ = ^ {{N2 - Nsf - 2Nf - \Ni\{N2 + N3)) , 
S- = ^ {{Nl - Ni) + \N,\{N2 - Ns)) . 

Since Eq. (14) with S_ — )■ and E'_ — t- implies that S"- — t- and since Ni — t- 0, 
it follows in the same way as in Bianchi VIIq that N2 — N3 ^ and A^| — N^ — )• 0. 
This means in particular that N2* = N^*- Eq. (14) with S4. — )• S+* 7^ and S^ — )■ 0, 
however, implies that 5*+ must not vanish asymptotically, but rather 

l\Ni\{N2 + N3) ^ {2 - q,)^+,. (46) 

A necessary condition for |A^i|(A^2 + A3) approaching a finite non-zero value is that the 
exponents of A^i, A'2 and N3 in Eq. (45) match: 

q^ — 4S+* + g* + 2S+* = <^ g* = T,^^. (47) 

Since we can write K according to Eq. (7) as 

K = ^ {{N2 - N^f + 2\Ni\{N2 + A^s)) , 

and hence conclude that A' — t- (2 — g*)S4_*, the Friedmann constraint implies 

1 = S^, + xl + yl + (2 - g,)S+,. (48) 

We get another equation from Eq. (18) with x' = 0, which yields 

x,{q,-2)-Foyl = 0. (49) 

We now solve the algebraic equations (47), (48), (49) for the unknowns S+*, x=k and y* 
where we replace q^: by Eq. (44). Under the condition that y* > and that also Eq. (19) 
with y' = must be satisfied, we find precisely one solution 

^ _ 4-2 _ /3 ci _ /3 Vc? + 2 



^* = \ n 2,1 ' y* = \ n 2,-1 - (50) 



2(c2 + l)' * V2cf + 1' '* \2 cl + l 
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Fig. 19: Deceleration scalar q (BVIII). 
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Fig. 20: Anisotropy variable S+ (BVIII). 
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Fig. 21: Decay of N^ (BVIII). 



Fig. 22: Evolution of N2/N3 (BVIII). 



This implies that 



c2 _2 



and hence this represents a decelerated epoch if and only if 

ci > \/2. 



As an additional check, we compute the exponent of A^^i in Eq. (45) in order to confirm 
that it is indeed negative 

'^* ~ ^^+* - 2 iT4- 

Hence, we claim that if ci > v2, generic Bianchi VIII solutions behave like Eqs. (43), 
(45) with A''2* = N^^, (46) and (50) asymptotically for r — >• 00. 

Numerical studies. For our numerical studies, we choose ci = 1.45 and C2 = 1.0 
and three different choices of initial data. This implies that 5]-(_* ~ 0.016, x* ~ 0.57, 
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Fig. 23: Decay ofY.^ (BVIII). 



Fig. 24: Constraint violations G using the 
modified evolution equations (BVIII). 



y* ^ 0.80 and g* « 0.016. Indeed, the variables S±, x, y, Ni and V approach these 
stationary values generically; see Figs. 19 and 20. We remark that we would get an 
even better agreement with the theoretical predictions if we would have let the code run 
beyond r = 400. In any case, we find that generic initially inflationary Bianchi VIII 
solutions have a graceful exit from inflation. In Fig. 21, we show the decay of A^i which 
is in agreement with Eq. (45); the bold line in Fig. 21 represents the expected decay. We 
demonstrate that A'^2 — -^3 — )• in Fig. 22. Similar to Bianchi VIIq, Fig. 23 shows that 
the decay of S_ is oscillatory and that the decay rate is smaller than in the Bianchi I case 
(as represented by the bold line in the figure). Just as a final check that our numerical 
code works reliably, we also plot the constraint violations in Fig. 24 for these runs. This 
plot demonstrates that our constraint damping scheme (see Section 3.2) works well for a 
long time. At some point, the numerical errors are dominated by the rapid growth of A'^2 
and N^ and subtle cancellations which need to take place in order to compute 5+ and 
S- , and hence the constraint violations start to grow. The time, however, for which the 
numerical solution is reliable, can be increased by decreasing the size of the numerical 
time step. 

Results. The situation for Bianchi VIII is similar to Bianchi VIIq as far as the behavior 
of the unbounded quantities A'^2 and N^ is concerned. The other variables T,±, x, y, Ni 
and "0 are bounded and, in fact, approach stationary values asymptotically with certain 
universal values which only depend on ci . The Bianchi VIIq future attractor discussed in 
the previous section is not an attractor for Bianchi VIII solutions. Its behavior is quite 
different; for example, Bianchi VIII solutions do not isotropize and become spatially 
flat. In comparison of this and Theorem 1.2 in [30], it turns out that the asymptotic 
properties of generic Bianchi VIII solutions are similar in the vacuum and the scalar 
field case (in contrast to the Bianchi VIIq case). Indeed, the same asymptotic behavior 
as in the vacuum case is obtained formally in the limit ci — t- oo. In any case, our main 
claim that generic initially inflationary Bianchi VIII scalar field models have a graceful 
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exit from inflation (if ci is chosen appropriately) has been confirmed. 

4.5 Beyond the monotonic case of the potential 

4.5.1 Monotonic potentials vs. potentials with a valley 

We recall the three main cases for the scalar field potential V{(p) given by Eq. (31) in 
Fig. 4 which mainly differ in the number of critical points of the function V. Notice 
that the critical points correspond to the zeros of the function F{(j)) defined in Eq. (12) 
which are given by 

, 1 ± Vl - Cfc2 

(Pl,2 = • 

Cl 

There are two (real) critical points — a local minimum and a local maximum as for the 
third plot in Fig. 4 — if and only if 1 — c^C2 > 0. We shall refer to this case as the 
"potential with a valley" . The potential has only one critical point — a saddle point as 
in the second plot in the figure — if and only if 1 — cfc2 = 0. We shall ignore this case 
in the following. Finally, there are no critical points — as for the first plot in the figure 
— if and only if 1 — cf C2 < 0; this is the "monotonic" case which we have considered 
exclusively so far. In this latter case, we have seen that the scalar field <j) is able to roll 
down the potential hill all the way to infinity. This is why our asymptotic analysis based 
on the limit i?i> — t- oo (or equivalently -0 = l/(f) — t- 0) discussed in the previous sections is 
relevant. On the other hand, if the potential has a valley, it is expected to depend on 
the initial data if the scalar field reaches infinity. If it does, then we expect the same 
asymptotics as discussed for the monotonic potential. If it does not and is rather "stuck" 
in the valley, it is expected that the future asymptotics of the solutions is different. 

Let us consider the potential with a valley. Since -F(0i) = at the minimum (j) = (pi 
(as well as at the maximum), the condition x' = (due to stationarity) in Eq. (18) implies 
that X = (if we exclude the case q = 2 for y > for the same reason as before) . Then 
Eq. (19) together with y' = yields that q = —1. The remaining evolution equations 
hence yield the stationary solution 

S± = 0, Ni,N2,N3 = 0, x = 0, y = l, (/) = h, 

with q = —1. The standard stability analysis reveals that this stationary point is hy- 
perbolic and in particular future stable. Indeed, it is just the de-Sitter solution where 
the non-dynamic scalar field represents a positive cosmological constant. The results in 
[26] suggest that solutions for which the scalar field is trapped in the valley approach 
this stationary point asymptotically, and hence infiation never stops and is de-Sitter like 
asymptotically. This is indeed confirmed numerically below. A similar analysis for the 
maximum (f) = (p2 of the potential yields that there exists no stable stationary point — 
in agreement with expectations. 

4.5.2 Numerical studies of the potential with a valley 

We choose the parameters ci = 1.45, C2 = 0.3 which corresponds to the case of two 
critical points, i.e., the "case with a valley". The two critical points represent the 
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locations of a local minimuni at (pi ~ 0.27 and a local maximum at 4>2 ~ 1.11. For 
this discussion, we restrict to the simplest Bianchi case — Bianchi I — and explore 
numerically whether solutions starting from inflationary initial data have graceful exits 
or not. 

In Fig. 25, we compute the projections of the orbits to the x-y-plane in the state 
space (similar to Fig. 9). We choose three different inflationary initial data, one where 
the initial value (po of the scalar field is 0.1 and hence satisfies (po < (pi (on the "downhill 
side" inside the valley), one where (pQ has the value 1.0 and hence satisfies 0i < ^o < 4*2 
(on the "uphill side" inside the valley) and one where (pQ has the value 2.0 and hence 
satisfies (p2 < (po (outside of the valley) . All three choices of initial data have the same 
positive initial value xq = 0.1 of x, and hence the scalar field increases initially (recall 
that X is proportional to (p). The outcome is in agreement with our expectations, see 
Fig. 25. The solution corresponding to the first choice of initial data (notice that the 
solutions corresponding to all three initial data start from the same point in the x-y-plane 
of the state space) has the property that the scalar field roles down the downhill side of 
the valley and hence x increases initially. It crosses the minimum of the potential with a 
positive (almost maximal velocity x) and then roles up the uphill side inside the valley. 
This has the consequence that x decreases and eventually becomes zero somewhere half 
way up this side of the valley. Then x becomes negative and the scalar field starts to 
role down towards the minimum again. After a while the whole process repeats but with 
a smaller velocity. Because q is negative during this whole process, our investigations 
in Section 4.2 show that the orbits must approach the circle x^ + y^ = 1 in Fig. 25 
and hence isotropize. Hence during each period of oscillation around the minimum, the 
solution gets closer to the stationary point associated with the minimum above, i.e., the 
de-Sitter solution. This solution is therefore never able to escape the infiationary regime 
and hence there is no graceful exit. A very similar dynamics is yielded for the second 
choice of initial data above. The only difference is that the scalar field roles up the uphill 
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side of the valley at the beginning. However, because the initial velocity x is not large 
enough, it is not able to escape from the valley and therefore has the same destiny as 
in the first case above: eventually the solution approaches the stationary point which 
represents the de-Sitter solution and there is therefore no graceful exit. In the third case, 
the scalar field is already outside of the value initially and initially moves away from the 
valley towards infinity. There is hence nothing which could stop the scalar field from 
reaching infinity and hence from having the same asymptotics as discussed in the case of 
the monotonic potential in previous sections. This is confirmed by comparing the third 
curve in Fig. 25 with the curves in Fig. 9. 

In Fig. 26 we study the case of three different sets of inflationary initial data, all 
with the same <j)o = 1-0 (i.e., on the uphill side inside of the valley). This time, however, 
we vary the initial value xq of x, i.e., the initial velocity of the scalar field. If the initial 
velocity is positive but small, the solution is trapped inside the valley as before and 
hence inflation never stops. If the initial velocity is positive and sufficiently large, we 
expect that the solution is able to escape the valley and hence reach infinity. This is 
confirmed in Fig. 26. In particular, we notice that then the future asymptotics is the 
same as in the case of the strictly monotonic potential: the solution approaches the same 
stationary solution as above, and, in particular, there is a graceful exit from inflation. 

5 Discussion 

In this paper we have investigated cosmological models with graceful exits from inflation. 
We have restricted to minimally coupled self-gravitating Bianchi A scalar field models 
with a particular form of the potential. We have chosen the potential specifically in 
order to violate known sufficient conditions for eternal infiation in the literature. Most 
of the studies of the graceful exit problem in the literature are restricted to the spatially 
homogeneous and isotropic case; see for example [25], where the graceful exit problem 
has been considered for a very similar class of potentials (in Eq. (31), one sets C2 = 
and the factor 0^ is replaced by (^" for a general n). We find that the presence of fu- 
ture attractors in the decelerated regime, if the parameter ci of the potential is chosen 
specifically and if the scalar field "escapes to infinity" (which is always the case if our 
potential is monotonic), guarantees that initially infiating models become decelerated 
after a finite evolution time. In particular, they approach one of these attractors. In the 
Bianchi cases I, II, VIq, these future attractors are stationary solutions of the dynam- 
ical system for the Hubble normalized variables. In the Bianchi cases VHq and VIII, 
however, some of the geometric variables are unbounded and hence the asymptotic be- 
havior is more difficult. Nevertheless, we are able to describe the asymptotic behavior 
of these solutions. Our investigations yield remarkable conclusions about the process 
of isotropization. While it is straightforward to prove that our models must isotropize 
and the spatial curvature must decay during infiation, it depends on the Bianchi case 
whether this process continues or stops when inflation is over. In fact, for the Bianchi 
cases II and VIq and VIII, this process stops and anisotropies and spatial curvature grow 
again after inflation. If the potential is not monotonic but rather has a "valley" then it 
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depends on the particular choice of the initial conditions whether a graceful exit from 
inflation occurs. 

As we have noticed before, we believe that some of the results in this paper generalize 
to more general classes of potentials; a particularly important example is the exponential 
potential. In fact, only little information about the potential is used in the analysis of 
the problem. This conjecture, however, needs to be studied in detail in future work. 
Also more work is required to understand the oscillatory behavior and the particular 
rates of decay in the Bianchi VIIq and VIII cases. 

Our techniques also allow us to determine the time between the begin of inflation 
and the graceful exit for our models. Physical arguments [24] lead to the conclusion that 
inflation should end after 75 — 100 efolds, which corresponds to the time r = 75 — 100 for 
the graceful exit. The actual time when inflation stops for our models certainly depends 
on the particular choice of parameters and initial data. However, our investigations 
show that it is not necessary to "fine-tune" the models in order to achieve a graceful exit 
around this time. Another interesting related question is about the accelerated expansion 
of the present universe due to dark energy which is suggested by current observations; 
see [22] for recent observational data and [24] for the theoretical background. The 
models, which we have presented in our paper, are not compatible with this since the 
decelerated epoch after the graceful exit lasts forever. However, as suggested by the 
results in [26] and confirmed by further numerical investigations, which we did, we find 
that if we add an arbitrarily small positive constant 14 to our potential in Eq. (31) 
(which is analogous to adding a positive cosmological constant to Einstein's equations), 
the corresponding solutions behave like the models presented in this paper for a long time 
of the evolution (in particular, there is a graceful exit from inflation). Then, however, 
at some late time which depends on the size of 14, these models suddenly deviate from 
this and eventually approach the de-Sitter solution. Hence, by this slight modiflcation, 
we are able to model the three main phases of the history of our universe: inflation, 
the intermediate decelerated epoch and the late time accelerated epoch of our present 
universe. 

Our investigations are mainly based on numerical methods. Due to a constraint 
violating unstable mode near the solutions of interest and the particular form of the 
subsidiary system, it has been crucial for our investigations to make sure that the Fried- 
mann constraint is satisfied as accurately as possible during the evolution. We have 
achieved this by introducing constraint damping terms to the evolution equations. 

In this whole paper, we have excluded the Bianchi IX case. The main reason is that it 
is not guaranteed that the Hubble normalized variables, which we use, are well-defined. 
In fact, the Bianchi IX case should be studied independently, for example, in terms of 
the variables introduced in [17]. 
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